clear; close all

qbit=1;showfig = 0;

Nsub=256;   % no. of subcarriers for rate computation
Ngrid = 1000;
L = 10e-1;

factor = 2.5;
x1 = 40e-1*factor;
y1 = 50e-1*ones(size(factor));
z1 = 90e-1*factor;
x2 = -30e-1*factor;
y2 = 50e-1*ones(size(factor));
z2 = 40e-1*factor;
TRX = zeros(6,length(factor));
TRX(1,:) = x1;
TRX(2,:) = y1;
TRX(3,:) = z1;
TRX(4,:) = x2;
TRX(5,:) = y2;
TRX(6,:) = z2;
load('TRX_std.mat','TRX','NTRX')
NTRX = 200;
dist = rms(TRX(1:3,:))*sqrt(3)+rms(TRX(4:6,:))*sqrt(3);

fc = 30e9;
c = 3e8;
BW_des =  [1e-8 0.2:0.2:3]*1e9;
Ptx = 10;%dBm

cal_rate = 0;
if cal_rate

Rates = zeros(3,length(Ptx),NTRX);
Ratessub = Rates;
RatesGS = Rates;
for n = 1:length(BW_des)
    BW = BW_des(n);
    nsub = max(Nsub,round(BW/1e9*200));
    parfor m = 1:NTRX
    Rates(:,n,m) = my_compute_rate(10,nsub,L,TRX(:,m),fc,BW,Ptx,showfig);
    Ratessub(1,n,m) = compute_rate_sub(10,nsub,L,TRX(:,m),fc,BW,Ptx,showfig,[2 0]);
    RatesGS(1,n,m) = compute_rate_GS(10,nsub,L,TRX(:,m),fc,BW,Ptx,showfig);
    end
    disp(n)
end
Rates = mean(Rates,3);
Ratessub = mean(Ratessub,3);
RatesGS = mean(RatesGS,3);
save('data/Rates_sum_rate_vs_bandwidth.mat','Rates','BW_des','Ratessub','RatesGS')

else 
    load('data/Rates_sum_rate_vs_bandwidth.mat')
end

BW_des = BW_des/1e9;
figure(1);

plot(BW_des,Rates(3,:),'--','LineWidth',2) % optimal
hold on
grid on
plot(BW_des,Rates(2,:),'b*-','LineWidth',2) % proposed
plot(BW_des,RatesGS(1,:),'rv-.','LineWidth',2) % proposed
plot(BW_des,Rates(1,:),'k^--','LineWidth',2); % classial 
plot(BW_des,Ratessub(1,:),'+--','LineWidth',2); % subarray 
xlabel('Bandwidth (GHz)')
ylabel('Average achievable rate (bit/s)')
% legend('Optimal beamforming','Proposed method','Proposed method, b = 3' ...
%     ,'Proposed method, b = 2','Proposed method, b = 1','Classical beamforming','Classical beamforming, b = 1' ...
%     ,'Location','southeast')
legend('Optimal beamforming','FZ-SPM','FZ-GSA','Classical beamforming',...
    'VSA-based beamforming',...
    'Location','southeast', 'FontSize', 14)
hold off
set(get(gca, 'Legend'), 'Box', 'off'); 
set(gca, 'FontSize', 16);
set(gcf, 'Position', [100 100 600 500]); % 调整当前图窗的大小
set(gca,'FontName','Times New Roman')
set(0, 'DefaultAxesFontName', 'Times New Roman');
saveas(gcf,'fig/Fig_sum_rate_vs_bandwidth.fig')
saveas(gcf,'fig/Fig_sum_rate_vs_bandwidth.png')

figure(2);
plot(BW_des,Rates(3,:),'--','LineWidth',1.5) % optimal
hold on
% plot(BW_des,Rates(2,:),'v-','LineWidth',1.5) % proposed
plot(BW_des,Rates(1,:),'k^--','LineWidth',1.5); % classial 
xlabel('Bandwidth (GHz)')
ylabel('Average rate per subcarrier (bit/s)')
set(gca,'FontName','Times New Roman')
% legend('Optimal beamforming','Proposed method','Proposed method, b = 3' ...
%     ,'Proposed method, b = 2','Proposed method, b = 1','Classical beamforming','Classical beamforming, b = 1' ...
%     ,'Location','southeast')
legend('Optimal beamforming','Classical beamforming',...
    'Location','southeast')
hold off
